﻿using System;

namespace DegreeDistance
{
    /// <summary>
    /// 经纬度坐标
    /// </summary>    

    public class Degree
    {
        public double X
        {
            get;
            set;
        }
        public double Y
        {
            get;
            set;
        }

        public override string ToString()
        {
            return string.Format("[x={0}\ty={1}]", X, Y);
        }
    }


    public class CoordDispose
    {
        private const double EarthRadius = 6378137.0d;//地球半径(米)

        /// <summary>
        /// 角度数转换为弧度公式
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        private static double Radians(double d)
        {
            return d * Math.PI / 180.0;
        }

        /// <summary>
        /// 弧度转换为角度数公式
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        private static double Degrees(double d)
        {
            return d * (180 / Math.PI);
        }

        /// <summary>
        /// 计算两个经纬度之间的直接距离
        /// </summary>
        public static double GetDistance(Degree degree1, Degree degree2)
        {
            double radLat1 = Radians(degree1.X);
            double radLat2 = Radians(degree2.X);
            double a = radLat1 - radLat2;
            double b = Radians(degree1.Y) - Radians(degree2.Y);

            double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) +
             Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2)));
            s = s * EarthRadius;
            s = Math.Round(s * 10000) / 10000;
            return s;
        }

        /// <summary>
        /// 计算两个经纬度之间的直接距离(google 算法)
        /// </summary>
        public static double GetDistanceGoogle(Degree degree1, Degree degree2)
        {
            double radLat1 = Radians(degree1.X);
            double radLng1 = Radians(degree1.Y);
            double radLat2 = Radians(degree2.X);
            double radLng2 = Radians(degree2.Y);

            double s = Math.Acos(Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Cos(radLng1 - radLng2) + Math.Sin(radLat1) * Math.Sin(radLat2));
            s = s * EarthRadius;
            s = Math.Round(s * 10000) / 10000;
            return s;
        }

        /// <summary>
        /// 以一个经纬度为中心计算出四个顶点
        /// </summary>
        /// <param name="degree1"></param>
        /// <param name="distance">半径(米)</param>
        /// <returns></returns>
        public static Degree[] GetDegreeCoordinates(Degree degree1, double distance)
        {
            double dlng = 2 * Math.Asin(Math.Sin(distance / (2 * EarthRadius)) / Math.Cos(degree1.X));
            dlng = Degrees(dlng);//一定转换成角度数  原PHP文章这个地方说的不清楚根本不正确 后来lz又查了很多资料终于搞定了

            double dlat = distance / EarthRadius;
            dlat = Degrees(dlat);//一定转换成角度数

            return new[]
                       {
                           //left-top
                           new Degree {X = Math.Round(degree1.X + dlat, 6), Y = Math.Round(degree1.Y - dlng, 6)},
                           //left-bottom
                           new Degree {X = Math.Round(degree1.X - dlat, 6), Y = Math.Round(degree1.Y - dlng, 6)},
                           //right-top
                           new Degree {X = Math.Round(degree1.X + dlat, 6), Y = Math.Round(degree1.Y + dlng, 6)},
                           //right-bottom
                           new Degree {X = Math.Round(degree1.X - dlat, 6), Y = Math.Round(degree1.Y + dlng, 6)}
                       };

        }
    }

    class Mercator2
    {
        //const double PI = 3.1415926535897932;//圆周率
        private const double PI = Math.PI;
        int IterativeTimes;     //反向转换程序中的迭代次数
        double IterativeValue;  //反向转换程序中的迭代初始值

        double EA;        //椭球体长半轴,米
        double EB;        //椭球体短半轴,米
        double EB0;        //基准纬度,弧度
        double EL0;        //原点经度,弧度

        public Mercator2()
        {
            IterativeTimes = 10;     //迭代次数为10
            IterativeValue = 0;      //迭代初始值
            EB0 = 0;
            EL0 = 0;
            //EA = 1;
            //EB = 1;

            //赤道半径为6378137米
            EA = 6378137.0;
            EB = 6356752.3142;
        }

        //角度到弧度的转换
        double DegreeToRad(double degree)
        {
            return PI * ((double)degree / (double)180);
        }

        //弧度到角度的转换
        double RadToDegree(double rad)
        {
            return (180 * rad) / PI;
        }

        //设定EA与EB
        void SetAB(double a, double b)
        {
            if (a <= 0 || b <= 0)
            {
                return;
            }
            EA = a;
            EB = b;
        }

        //设定EB0
        void SetB0(double b0)
        {
            if (b0 < -PI / 2 || b0 > PI / 2)
            {
                return;
            }
            EB0 = b0;
        }
        //设定EL0
        void SetL0(double l0)
        {
            if (l0 < -PI || l0 > PI)
            {
                return;
            }
            EL0 = l0;
        }

        //#############################################################################################
        //投影正向转换程序
        //double L: 经度,弧度
        //double B: 维度,弧度        
        //ref double X:   横向直角坐标
        //ref double Y:   纵向直角坐标
        //#############################################################################################
        public bool ToProj(double L, double B, ref double X, ref double Y)
        {
            L = DegreeToRad(L);
            B = DegreeToRad(B);
            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;
            double E = Math.Exp(1);
            if (L < -PI || L > PI || B < -PI / 2 || B > PI / 2)
            {
                return false;
            }

            f = (EA - EB) / EA;

            dtemp = 1 - (EB / EA) * (EB / EA);
            if (dtemp < 0)
            {
                return false;
            }
            e = Math.Sqrt(dtemp);

            dtemp = (EA / EB) * (EA / EB) - 1;
            if (dtemp < 0)
            {
                return false;
            }
            e_ = Math.Sqrt(dtemp);

            NB0 = ((EA * EA) / EB) / Math.Sqrt(1 + e_ * e_ * Math.Cos(EB0) * Math.Cos(EB0));

            K = NB0 * Math.Cos(EB0);
            X = K * (L - EL0);
            Y = K * Math.Log(Math.Tan(PI / 4 + B / 2) * Math.Pow((1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B)), e / 2));

            return true;
        }

        ///#############################################################################################
        //投影反向转换程序
        //double X: 纵向直角坐标
        //double Y: 横向直角坐标
        //ref double B: 维度,弧度
        //ref double L: 经度,弧度
        //#############################################################################################
        public bool FromProj(double pmtX, double pmtY, ref double pmtL, ref double pmtB)
        {
            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;
            double E = Math.Exp(1);

            if (EA <= 0 || EB <= 0)
            {
                return false;
            }

            f = (EA - EB) / EA;

            dtemp = 1 - (EB / EA) * (EB / EA);
            if (dtemp < 0)
            {
                return false;
            }
            e = Math.Sqrt(dtemp);

            dtemp = (EA / EB) * (EA / EB) - 1;
            if (dtemp < 0)
            {
                return false;
            }
            e_ = Math.Sqrt(dtemp);

            NB0 = ((EA * EA) / EB) / Math.Sqrt(1 + e_ * e_ * Math.Cos(EB0) * Math.Cos(EB0));

            K = NB0 * Math.Cos(EB0);

            pmtL = pmtX / K + EL0;
            pmtB = IterativeValue;

            for (int i = 0; i < IterativeTimes; i++)
            {
                pmtB = PI / 2 - 2 * Math.Atan(Math.Pow(E, (-pmtY / K)) * Math.Pow(E, (e / 2) * Math.Log((1 - e * Math.Sin(pmtB)) / (1 + e * Math.Sin(pmtB)))));
            }

            pmtL = RadToDegree(pmtL);
            pmtB = RadToDegree(pmtB);

            return true;
        }
    }
}
